 use "PNAS_AassveEtAl2024.dta",clear

**# Distribution and Mean/SD of Dependent Variables
* A_successfam is the DV of Experiment 1 & B_successfam is the DV of Experiment 2
 
  twoway (histogram A_successfam [fweight = weight], width(1) color(purple%30)) ///
       (histogram B_successfam [fweight = weight],  width(1) color(red%30)), ///
	   legend(order(1 "Experiment1" 2 "Experiment2" )) [fweight = weight] 

 tabstat A_successfam, by(country) stat(mean sd N) col(statistics) 

 svyset, clear
 svyset pid [pweight = weight], vce(linear)    

 svy: mean A_successfam 
 estat sd
 svy: mean B_successfam 
 estat sd
 
 svy: mean A_successfam, over(COUNTRY)
 estat sd
 svy: mean B_successfam, over(COUNTRY)
 estat sd


	   
**# Correlation betwen Factors/Demographic variables

 pwcorr union-WLB if experiment==1, sig
 pwcorr union-finsupport if experiment==2, sig

 pwcorr union-finsupport female edu_self married single catholic haschild working_self hhinc_adjusted if experiment==1, sig 
 pwcorr union-WLB female edu_self married single catholic haschild working_self hhinc_adjusted if experiment==2, sig 
 
 ************************************************************************
 ***                                                                  ***  
 ***             Main Effect Analysis - Random Intercept              ***
 ***                                                                  ***
 ************************************************************************
 
**# pooled 8 country: analysis

eststo clear
 eststo est1A: mixed A_successfam i.(union-WLB)  ib2.(income respect famcom extfamcom genderrole)  ib4.WLB i.run || country:  || ResponseId: [pweight = weight], nolog mle vce(robust)
 eststo est4A: mixed B_successfam i.(union-finsupport) ib2.(income respect famcom extfamcom genderrole finsupport childedu)  ib4.WLB i.run || country: || ResponseId: [pweight = weight], nolog mle vce(robust)


**# pooled 8 country: coefplots
 
 //set scheme s2color
 
 set scheme cleanplots
 coefplot (est1A, label(SuccessFam) mcolor("0 0 255") ciopts(color("0 0 255")) msymbol(O) )   ///
  || , drop(_cons *.run) xline(0) norecycle /// 
        headings(2.union = `""{bf: Union Status}" "{it: (Ref: married)}""'     ///
          2.fertility = `""{bf: Fertility}" "{it: (Ref: 1 child)}""'     ///
		  1.income = `""{bf: HH Income}" "{it: (Ref: average)}""'    ///
		  1.respect = `""{bf: Community respect}" "{it: (Ref: Not respected)}""'    ///
		  1.famcom = `""{bf: Fam. communication}" "{it: (Ref: Not well)}""'     ///
		  1.extfamcom = `""{bf: Ext. fam. contact}" "{it: (Ref: Not frequently)}""'   ///
	      1.genderrole = `" "{bf: Gender roles}" "{it: (Ref: woman }" "{it: double burden)}" "'    ///
		  1.WLB = `""{bf: Work-family conflict}" "{it: (Ref: Both conflicted)}""', gap(0.5) labsize(vsmall) ) /// 
		xlabel(,labsize(vsmall)) ylabel(,labsize(vsmall)) msize(vsmall)  legend(row(1) size(vsmall) region(lstyle(none))) ///
		grid(between glpattern(dot) glwidth(*0.5) glcolor(gray)) ///
		graphregion(margin(zero)) plotregion(margin(zero)) legend(off)
				    
  coefplot (est4A, label(SuccessFam) mcolor("0 0 255") ciopts(color("0 0 255")) msymbol(O) ) ///
  || , drop(_cons *.run) xline(0) norecycle /// 
        headings(2.union = `""{bf: Union Status}" "{it: (Ref: married)}""'     ///
          2.fertility = `""{bf: Fertility}" "{it: (Ref: 1 child)}""'     ///
		  1.income = `""{bf: HH Income}" "{it: (Ref: average)}""'    ///
		  1.respect = `""{bf: Community respect}" "{it: (Ref: Not respected)}""'    ///
		  1.famcom = `""{bf: Fam. communication}" "{it: (Ref: Not well)}""'     ///
		  1.extfamcom = `""{bf: Ext. fam. contact}" "{it: (Ref: Not frequently)}""'   ///
	      1.genderrole = `" "{bf: Gender roles}" "{it: (Ref: woman }" "{it: double burden)}" "'    ///
		  1.WLB = `""{bf: Work-family conflict}" "{it: (Ref: Both conflicted)}""'  ///
		  1.childedu = `""{bf: Child Educ.}" "{it: (Ref: Bachelors)}""'  ///
		  1.finsupport = `""{bf: Child fin. support}" "{it: (Ref: Not saving)}""' , gap(0.5)  labsize(vsmall) ) ///
 		xlabel(-1(0.5)0.5,labsize(vsmall)) ylabel(,labsize(vsmall)) msize(vsmall)  legend(row(1) size(vsmall) region(lstyle(none))) ///
		grid(between glpattern(dot) glwidth(*0.5) glcolor(gray)) ///
		graphregion(margin(zero)) plotregion(margin(zero)) legend(off)
		
 graph display, ysize(10) 
 
**# By country analysis and plots: Experiment 1

sort country 
eststo clear
by country: eststo: mixed A_successfam  i.(union-WLB) ib2.(income respect famcom extfamcom genderrole) ib4.(WLB) i.run  || ResponseId: [pweight = weight],  nolog vce(robust)


 esttab using result/bycountry_weighted.tex, replace se la nodep obslast starlevels(* 0.05 ** 0.01) ///
         transform(ln*: exp(@) exp(@)) ///
		 eqlabels("" "var(person)" "var(residual)", none) ///
		 varlabels(,elist(weight:_cons "{break}{hline @width}" "{break}{hline @width}")) ///
		 title(experiment A Successful Family Mixed-effect Model) 
	
	
 coefplot (est1, label(China) mcolor("255 0 0") ciopts(color("255 0 0")) msymbol(O) ) ///
		  (est2, label(Italy) mcolor("0 0 255") ciopts(color("0 0 255")) msymbol(O) ) ///
		  (est3, label(Japan) mcolor("255 135 255") ciopts(color("255 135 255")) msymbol(O) ) ///
		  (est4, label(Korea) mcolor("255 230 0") ciopts(color("255 230 0")) msymbol(O) ) ///
		  (est5, label(Norway) mcolor("100 180 255") ciopts(color("100 180 255")) msymbol(O) ) ///
		  (est6, label(Singapore) mcolor("255 140 0") ciopts(color("255 140 0")) msymbol(O) ) ///
		  (est7, label(Spain) mcolor("0 150 0") ciopts(color("0 150 0")) msymbol(O) ) ///
		  (est8, label(USA) mcolor("165 115 255") ciopts(color("165 115 255")) msymbol(O) ) ///
 || , drop(_cons *.run) xline(0)  ///
 order(*union *fertility *income *respect *famcom *extfamcom *genderrole *WLB )  ///
 xlabel(,labsize(vsmall)) ylabel(,labsize(vsmall)) msize(vsmall)  ///
 legend( position(6) row(2) size(vsmall) region(lstyle(none))) ///
 headings(2.union = `""{bf: Union Status}" "{it: (Ref: married)}""'     ///
          2.fertility = `""{bf: Fertility}" "{it: (Ref: 1 child)}""'     ///
		  1.income = `""{bf: HH Income}" "{it: (Ref: average)}""'    ///
		  1.respect = `""{bf: Community respect}" "{it: (Ref: Not respected)}""'    ///
		  1.famcom = `""{bf: Fam. communication}" "{it: (Ref: Not well)}""'     ///
		  1.extfamcom = `""{bf: Ext. fam. contact}" "{it: (Ref: Not frequently)}""'   ///
	      1.genderrole = `" "{bf: Gender roles}" "{it: (Ref: woman }" "{it: double burden)}" "'    ///
		  1.WLB = `""{bf: Work-family conflict}" "{it: (Ref: Both conflicted)}""'   , gap(0.5)  labsize(vsmall) ) ///
 graphregion(margin(0 +0.5 0 0)) plotregion(margin(zero)) 
 graph display, ysize(12)
 
**# By country analysis and plots: Experiment 2

sort country 
eststo clear
by country: eststo: mixed B_successfam  i.(union-finsupport) ib2.(income respect famcom extfamcom genderrole childedu finsupport) ib4.(WLB) i.run  || ResponseId: [pweight = weight],  nolog vce(robust)


		 
 coefplot (est1, label(China) mcolor("255 0 0") ciopts(color("255 0 0")) msymbol(O) ) ///
		  (est2, label(Italy) mcolor("0 0 255") ciopts(color("0 0 255")) msymbol(O) ) ///
		  (est3, label(Japan) mcolor("255 135 255") ciopts(color("255 135 255")) msymbol(O) ) ///
		  (est4, label(Korea) mcolor("255 230 0") ciopts(color("255 230 0")) msymbol(O) ) ///
		  (est5, label(Norway) mcolor("100 180 255") ciopts(color("100 180 255")) msymbol(O) ) ///
		  (est6, label(Singapore) mcolor("255 140 0") ciopts(color("255 140 0")) msymbol(O) ) ///
		  (est7, label(Spain) mcolor("0 150 0") ciopts(color("0 150 0")) msymbol(O) ) ///
		  (est8, label(USA) mcolor("165 115 255") ciopts(color("165 115 255")) msymbol(O) ) ///
 || , drop(_cons *.run) xline(0)  ///
 order(*union *fertility *income *respect *famcom *extfamcom *genderrole *WLB )  ///
 xlabel(,labsize(vsmall)) ylabel(,labsize(vsmall)) msize(vsmall)  ///
 legend( position(6) row(2) size(vsmall) region(lstyle(none))) ///
 headings(2.union = `""{bf: Union Status}" "{it: (Ref: married)}""'     ///
          2.fertility = `""{bf: Fertility}" "{it: (Ref: 1 child)}""'     ///
		  1.income = `""{bf: HH Income}" "{it: (Ref: average)}""'    ///
		  1.respect = `""{bf: Community respect}" "{it: (Ref: Not respected)}""'    ///
		  1.famcom = `""{bf: Fam. communication}" "{it: (Ref: Not well)}""'     ///
		  1.extfamcom = `""{bf: Ext. fam. contact}" "{it: (Ref: Not frequently)}""'   ///
	      1.genderrole = `" "{bf: Gender roles}" "{it: (Ref: woman }" "{it: double burden)}" "'    ///
		  1.WLB = `""{bf: Work-family conflict}" "{it: (Ref: Both conflicted)}""'  ///
		  1.childedu = `""{bf: Child Educ.}" "{it: (Ref: Bachelors)}""'  ///
		  1.finsupport = `""{bf: Child fin. support}" "{it: (Ref: Not saving)}""' , gap(0.5)  labsize(vsmall) ) ///
  graphregion(margin(0 +0.5 0 0)) plotregion(margin(zero))  
 graph display, ysize(15) 
 
**# Robustness check
* Experiment 1
eststo clear
 eststo: mixed A_successfam i.(union-WLB)  ib2.(income respect famcom extfamcom genderrole)  ib4.WLB i.run || country:  || ResponseId: [pweight = weight], nolog mle vce(robust)
  eststo: mixed A_successfam i.(union-WLB)  ib2.(income respect famcom extfamcom genderrole)  ib4.WLB i.run || country:  || ResponseId: [pweight = weight] if age!=., nolog mle vce(robust)
 eststo: mixed A_successfam i.(union-WLB)  ib2.(income respect famcom extfamcom genderrole)  ib4.WLB i.run || country: i.(union-WLB) || ResponseId: [pweight = weight], nolog mle vce(robust)
 xtset COUNTRY
 eststo: reghdfe A_successfam i.(union-WLB)  ib2.(income respect famcom extfamcom genderrole)  ib4.WLB i.run [pweight = weight], a(COUNTRY) vce(robust)
 eststo : xtreg A_successfam i.(union-WLB)  ib2.(income respect famcom extfamcom genderrole)  ib4.WLB i.run , fe vce(robust)

* Experiment 2
 eststo est5: mixed B_successfam i.(union-finsupport) ib2.(income respect famcom extfamcom genderrole finsupport childedu)  ib4.WLB i.run || country: || ResponseId: [pweight = weight], nolog mle vce(robust)
 eststo est6: mixed B_successfam i.(union-finsupport) ib2.(income respect famcom extfamcom genderrole finsupport childedu)  ib4.WLB i.run || country: i.(union-finsupport) || ResponseId: [pweight = weight], nolog mle vce(robust)
 eststo est7: reghdfe B_successfam i.(union-finsupport) ib2.(income respect famcom extfamcom genderrole finsupport childedu)  ib4.WLB i.run [pweight = weight], a(COUNTRY) vce(robust)
 eststo est8: xtreg B_successfam i.(union-finsupport) ib2.(income respect famcom extfamcom genderrole finsupport childedu)  ib4.WLB i.run, fe vce(robust)

**# Descrpitve statistics - wide data 
		
use  "PNAS_AassveEtAl2024_wide.dta", clear
**#  unweighted

 tabstat age female highschool higheredu_self single married cohabiting nchild* hhinc_adj, by(country) stat(mean sd N)
 
**#  weighted

 svyset, clear
 svyset ResponseId [pweight = weight], vce(linear)    

 svy: mean age , cformat(%9.4f)
 estat sd 
 svy: mean female , cformat(%9.4f)
 estat sd
 svy: mean highschool , cformat(%9.4f)
 estat sd
 svy: mean higheredu_self , cformat(%9.4f)
 estat sd
 svy: mean single , cformat(%9.4f)
 estat sd
 svy: mean married , cformat(%9.4f)
 estat sd
 svy: mean cohabiting , cformat(%9.4f)
 estat sd
 svy: mean nchild_wzero , cformat(%9.4f)
 estat sd
 svy: mean nchildren , cformat(%9.4f)
 estat sd
 svy: mean nchilddum1 , cformat(%9.4f)
 estat sd
 svy: mean nchilddum2 , cformat(%9.4f)
 estat sd
 svy: mean nchilddum3 , cformat(%9.4f)
 estat sd
 svy: mean nchilddum4 , cformat(%9.4f)
 estat sd
 svy: mean hhinc_adj , cformat(%9.4f)
 estat sd

 svy: mean age country, col percent 
 estat sd 

 svy: mean female country, col percent 
 estat sd 

 svy: mean highschool ,over(COUNTRY) cformat(%9.2f)
 estat sd
 svy: mean higheredu_self ,over(COUNTRY) cformat(%9.2f)
 estat sd

 svy: mean single ,over(COUNTRY) cformat(%9.2f)
 estat sd
 svy: mean married ,over(COUNTRY) cformat(%9.2f)
 estat sd
 svy: mean cohabiting ,over(COUNTRY) cformat(%9.2f)
 estat sd
 
 svy: mean nchild_wzero ,over(COUNTRY) cformat(%9.2f)
 estat sd
 svy: mean nchildren ,over(COUNTRY) cformat(%9.2f)
 estat sd

 svy: mean  nchilddum1 ,over(COUNTRY) cformat(%9.2f)
 estat sd
 svy: mean  nchilddum2 ,over(COUNTRY) cformat(%9.2f)
 estat sd
 svy: mean  nchilddum3 ,over(COUNTRY) cformat(%9.2f)
 estat sd
 svy: mean  nchilddum4 ,over(COUNTRY) cformat(%9.2f)
 estat sd

 svy: mean hhinc_adj ,over(COUNTRY) cformat(%9.2f)
 estat sd
 
 





